Skip to content
This repository has been archived by the owner on Dec 22, 2021. It is now read-only.

Stabilize derivatives in fatiando.gravmag #196

Merged
merged 15 commits into from May 26, 2015
Merged

Stabilize derivatives in fatiando.gravmag #196

merged 15 commits into from May 26, 2015

Conversation

leouieda
Copy link
Member

As pointed out by @mycarta in #194, the FFT-based derivative calculations in fatiando.gravmag are unstable (see the figure below for a comparison with the MATLAB derivatives).

deriv

This PR tries to stabilize the derivatives by padding the input data and using different algorithms.

Fixes #194 and #167

Methods

The input data grids for FFT derivatives are padded using numpy.pad with the values in the borders of the grid. The padding is done to the nearest power of two. This solves the instabilities in the z derivative (specially in the edges). The horizontal derivatives are a bit better but still suffer from the "striping" seen above.

Horizontal derivatives can be stabilized by using a central differences scheme instead of the FFT. This is what matlab uses to get the stable derivatives seen in the figure. The derivx and derivy functions will have options to use the finite differences instead of FFT.

Checklist:

  • Make tests for new code
  • Create/update docstrings
  • Include relevant equations and citations in docstrings
  • Code follows PEP8 style conventions
  • Code and docs have been spellchecked
  • Include new dependencies in docs, requirements.txt, README, and .travis.yml
  • Documentation builds properly
  • All tests pass
  • Check cookbook recipes
  • Can be merged
  • Changelog entry (leave for last)
  • Firt-time contributor? Add yourself to doc/contributors.rst (leave for last)
  • Update I/O functions to read data with shape = (nx, ny)
  • Warn on the mailing list about the changes to gridder.regular

@leouieda
Copy link
Member Author

I found an error in the way fatiando.gridder.regular works. In Fatiando, we assume that x is North-South and y is East-West. So if data is on a matrix, it would make sense that lines are x and columns are y. In the gridder.regular function we were assuming the opposite. The functions asked for a shape = (ny, nx). The correct form should be shape = (nx, ny).

This might impact some functions but it will be restricted mostly to the FFT derivatives so I'll make the fix in this branch.

Padding helps reduce edge effects of the derivatives. Use the edge
values to pad the array to the nearest power of two (almost, actually).
Was using ny, nx = shape when in fact x is north-south so
shape = nx, ny
gridder.regular is generating grids the wrong way. The shape should be
(nx, ny), not (ny, nx). Thus, the derivatives are in the wrong direction
and the tests fail.
The Euler deconvolution tests were using FFT derivatives and failing
because of them. I exchanged for derivatives approximated by finite
differences in the forward modeling. This was the Euler tests don't
depend of the performance of the FFT derivs.
This is how it should be. It breaks the tests for gravmag.sphere which
were hardcoded numbers as doctests. Now the numbers are out of order so
the test fail. Need a better test for this.
These tests weren't helping much. They compared the output to previous
output. There are better tests comparing numpy and Cython
implementations in the unit tests.
Filled the missing arguments in the docstrings and make explicit using
fft derivs in the tests.
@leouieda
Copy link
Member Author

@mycarta here are the update results for this branch:

Old results

New results

index

Generated using this notebook.

Notice that the x and y derivatives are exchanged. There is a lot of confusion regarding this because we use x -> North-South in Fatiando and that is not always the case with other software. It's good to check that the results have the proper orientation that you'd expect.

For the horizontal derivatives, the finite difference solutions are much
more stable than the FFT ones. They underestimate the derivative a bit
where the gradient is high but otherwise behave in a more consistent
manner and are safer to use.
@mycarta
Copy link

mycarta commented Apr 30, 2015

Hi Leo

Thanks for the follow-up on derivatives.

Could you please take the coordinates off from the images, and the notebook out altogether or use
mock coordinates in it?

Those coordinates are from my unpublished thesis and not supposed to be out
there, they were for your eyes only.

My apologies I should have said it clearly.

In my tutorial I'd be using mock coordinates.

Thanks

Matteo

On Thu, Apr 30, 2015 at 9:30 AM, Leonardo Uieda notifications@github.com
wrote:

@mycarta https://github.com/mycarta here are the update results for
this branch:
Old results

https://cloud.githubusercontent.com/assets/290082/7272735/9a72e110-e8c3-11e4-907f-5876ded6bd8f.png
New results

[image: index]
https://cloud.githubusercontent.com/assets/290082/7416251/6ed571b8-ef34-11e4-801a-943adc275ea6.png

Generated using this notebook
http://nbviewer.ipython.org/gist/leouieda/ac346557bd28f77861b8.

Notice that the x and y derivatives are exchanged. There is a lot of
confusion regarding this because we use x -> North-South in Fatiando and
that is not always the case with other software. It's good to check that
the results have the proper orientation that you'd expect.


Reply to this email directly or view it on GitHub
#196 (comment).

@leouieda
Copy link
Member Author

I'm so sorry about that! I should have noticed that your plots have no coordinates. It's fixed now. Also removed from the notebook (which doesn't include any data).

@mycarta
Copy link

mycarta commented Apr 30, 2015

If I'd mentioned it to begin with.....
All's well that ends well
:)
Thanks Leo

On Thu, Apr 30, 2015 at 10:54 AM, Leonardo Uieda notifications@github.com
wrote:

I'm so sorry about that! I should have noticed that your plots have no
coordinates. It's fixed now. Also removed from the notebook (which doesn't
include any data).


Reply to this email directly or view it on GitHub
#196 (comment).

leouieda added a commit that referenced this pull request May 26, 2015
Stabilize derivatives in fatiando.gravmag
@leouieda leouieda merged commit 274007e into master May 26, 2015
@eusoubrasileiro
Copy link
Contributor

@leouieda awesome pictures of smootheness

mtb-za added a commit to mtb-za/fatiando that referenced this pull request May 27, 2015
Need to pull in some changes from upstream, particularly
fatiando#196 which has
made changes to `fatiando.gridder`, which this needs.
@leouieda leouieda deleted the fftderiv branch May 27, 2015 17:06
mtb-za added a commit to mtb-za/fatiando that referenced this pull request May 28, 2015
mtb-za added a commit to mtb-za/fatiando that referenced this pull request May 28, 2015
@leouieda leouieda mentioned this pull request Jul 6, 2015
11 tasks
Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Issues with gravmag fourier derivatives
3 participants